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Abstract. We study the spatial circular restricted prob- 
lem of three bodies in the light of Nekhoroshev theory of 
stability over large time intervals. We consider in partic- 
ular the Sun- Jupiter model and the Trojan asteroids in 
the neighborhood of the Lagrangian point L4. We find a 
region of effective stability around the point L4 such that 
if the initial point of an orbit is inside this region the orbit 
is confined in a slightly larger neighborhood of the equi- 
librium (in phase space) for a very long time interval. By 
combining analytical methods and numerical approxima- 
. tions we are able to prove that stability over the age of the 
<^> ' universe is guaranteed on a realistic region, big enough to 
O^l . include one real asteroid. By comparing this result with 
' the one obtained for the planar problem we sec that the 
regions of stability in the two cases are of the same mag- 
nitude. 

Key words: minor planets, asteroids - celestial mechan- 
ics - instabilities 
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1. Introduction 

The study of a Hamiltonian system in the neighborhood of 
an elliptic equilibrium point is of interest in many fields of 
mathematical physics and astronomy. Let us consider an 
analytic Hamiltonian H with n degrees of freedom, hav- 
ing an elliptic equilibrium point. Rigorous results prov- 
ing the existence of orbits which do not leave a neigh- 
borhood of the equilibrium can be given in the frame of 
KAM theory, under generic conditions of non-resonance 
and non-degeneracy. KAM guarantees the existence of 
many ri-dimensional invariant tori around the equilibrium 
point. However, such invariant tori do not fill an open re- 
gion, i.e. the possibility of the so-called Arnold diffusion 
cannot be excluded, except for the two dimensional case. 

An alternative approach is to look for results which are 
valid over a finite time interval, but give an effective bound 
on the Arnold diffusion. This goal can be achieved by con- 
structing the normal form of the Hamiltonian around the 
elliptic equilibrium point. Normal forms are a standard 
tool in Celestial Mechanics for studying the dynamics in 
the neighborhood of an elliptic equilibrium point. Usually 
these normal forms are obtained as divergent series but 
their truncation makes them useful. Roughly speaking one 
shows that the system admits a number of approximate in- 
tegrals, whose time variation can be controlled to be small 
for an extremely long time. In these cases we have effec- 
tive stability, i.e. even when an orbit is not stable, the time 
needed for it to leave the neighborhood of the equilibrium 
is larger than the expected lifetime of the physical system 
studied. This is the basis to derive the classical Nekhoro- 



realistic models of the solar system, and the statistical 



shev's estimates (Nekhoroshev |l977j) . Different proofs of 
t he N ekhoroshev theorem wer e giv en by Benettin et al. 
(|1985|), Benettin & Gallay otti (|l986|) , Giorgilli & Zehnder 
(|l992|) and Poschel (|1993| ). 

A scientific field where the Nekhoroshev theory has 
been applied is the problem of the stability of the Trojan 
asteroids. In recent years this problem has been investi- 
gated by a number of researchers, both numerically and 
analytically. The numerical investigations deal with the 
evolution in time of a sample of orbits, in sophisticated 



1993 



1994 



Levison et al. 



study of these orbit s (Mi lani 
19971 Tsiganis et al. |2000D . 

In analytical studies simpler models of the system have 
been used such as the two dimensional (2D) planar, and 
the three dimensional (3D) spatial restricted three body 
problem. According to Nekhoroshev theory, one has to 
estimate the rate of diffusion around the elliptic equilib- 
rium Lagrangian point L4. Because of the symmetries of 
the system the same study is valid for the L5 point. The 
problem has been previously investigated by Giorgilli et al. 
(1 19891 ), Celletti & Giorgilli ( |l99l| ) and Giorgilli & Skokos 
(1997) (hereafter Paper I). Also analytical stability results 
for particular orbits were provided by Celletti & Ferrara 
( 1996| ) who studied the orbit of the asteroid Ceres and by 
Jorba & Villanueva (1998) who studied a stable periodic 
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orbit around the equilibrium point L 5 . 

The estimation of the region of effective stability by 
Giorgilli et al. ( |l989| ) and Celletti & Giorgilli ( |l99l| ) was 
realistic but the region where the real asteroids were ac- 
tually found was larger by a factor 300 (in the best case) 
to 3,000, compared to the theoretical stability region. The 
theoretical estimation was significantly improved in Pa- 
per I, since the region found in the planar restricted three 
body problem was big enough to include 4 real asteroids, 
while most of them failed to be inside this region by a fac- 
tor only 10. This result underlined the fact that Nekhoro- 
shev theory can give meaningful estimates, applicable to 
real systems. On the other hand Arnold diffusion, which, 
as already mentioned, can drive some orbits with initial 
conditions near the equilibrium point L4 to regions of the 
phase space far away from it, appears only in the spatial 
case and not in the planar one. 

So in the present paper we study the spatial case, fol- 
lowing a similar procedure to the one used in Paper I, in 
order to find out to what degree the presence of Arnold 
diffusion changes the effective stability region. Apart from 
the fact that we consider the three dimensional case and 
not its restriction on the plane, some minor changes of the 
scheme used in Paper I improve slightly the estimation of 
the size of the stability region. In particular the expan- 
sion of the Hamiltonian of the system in a power series 
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suitable for the application of the normal form scheme is 
computed with higher accuracy than before. Also a more 
accurate calculation of the time needed for an orbit to 
leave a particular region of the phase space around the 
point L4 (escape time) is provided. We were also able to 
compute the normal form to higher orders than in previ- 
ous studies, both in the 3D and the 2D case. In particular 
in the 3D case we have numerically computed the normal 
form up to order 29, which is a hard task, since one has to 
manipulate functions with a much larger number of coef- 
ficients, than in the case of order 22 (Celletti & Giorgilli 



1991). In the 2D case we computed the normal form up 



to order 49 (instead of order 34 in Paper I). 

The paper is organized as follows. In Sect. 2 we present 
the Hamiltonian we use, referring to all the canonical 
transformations needed to bring it in a form suitable for 
the application of the normal form scheme. Also in this 
section we sketch the procedure of computing the normal 
form. The main results of the paper, concerning the esti- 
mation of the size of the regions of effective stability both 
in the spatial (3D) and the planar (2D) case, are presented 
in Sect. 3. The application of the above results to real as- 
teroids is also included in Sect. 3. In Sect. 4 we discuss 
some issues concerning the effectiveness of our scheme. 
Finally in Sect. 5 we summarize our results. 

2. The Hamiltonian and the normal form of the 
system 

The spatial restricted problem of three bodies, in partic- 
ular for the Sun (S), Jupiter (J) and asteroid (A) system 
can be described as follows: we study the motion of an 
asteroid A of infinitesimal mass, orbiting in the gravita- 
tional field of two primaries S and J with masses equal to 
1 — fi and fi respectively, which are assumed to revolve in 
circular orbits around their common center of mass. 

We introduce a uniformly rotating frame (O, q±, q 2 , 
q 3 ) so that its origin is located at the center of mass of the 
Sun- Jupiter system, with the Sun always at the point (fj,, 
0, 0) and Jupiter at the point (1 — fj,, 0, 0). The physical 
units are chosen so that the distance between Jupiter and 
the Sun is 1, fj, = 9.5387536 • 10~ 4 and the angular velocity 
of Jupiter is 1. The time unit is (2tt)^ 1 Tj, where Tj is the 
period of the circular motion of Jupiter around the Sun. 
So the age of the universe is about 10 10 time units. The 
Hamiltonian of the system is : 



q2Pi - qiP2 



1-/1 



V (qi - v) 2 + + il 



■92 + 93 



(1) 



VWi + 1 - m) 2 

The coordinates of the Lagrangian point L4 are: q\ 

M - 5i 12 = IT, 93 = 0, pi = p 2 = M ~ 5, P3 = 0. 

In order to bring the Hamiltonian to a form suitable for 
the application of the normal form scheme we perform a 



sequence of transformations. The first step is to introduce 
a uniformly rotating frame with its origin on the Sun (S) 
using the generating function 

W 3 = -(Qi + fi)pi - Q2P2 - Q3P3 + VQ2 , 

where Qi, Q2, Q3, Pi, P2, P3 are the heliocentric coordi- 
nates. 

It is known that the projection on the plane of Jupiter's 
orbit, of the stability region around L4, is a banana- 
shaped region which lies close to the circle with center the 
Sun and radius equal to the Sun-Jupiter distance. Since 
the plane of Jupiter's orbit is a symmetry plane for the 
system, a good choice for describing this region are the 
cylindrical coordinates P, 0, Z, which are introduced by 
the generating function 
W 3 = -P{Pi cos® + P 2 sinS) - ZP 3 . 
In this system of coordinates L4 is located at P — 1, 
9 = f , z = 0, vp = 0, Pe = l,Pz = 0. 

We move the origin of the coordinate system to the 
point L4 using the generating function 

2irp y 



W 2 = Px (P-l) + ( P y + l)G 

Then the Hamiltonian becomes: 



Pz 



H 



1 



-fj,(x- 



(x + 
1) cos 



I) 2 



\-p: 

2tt 

T 



Py 



y/(x+l) 2 + Z 2 



(2) 



yftx + l) 2 + z 2 + 1 + 2{x + 1) cos (y + 2*) 

We expand the above Hamiltonian in Taylor series 
around the point L4 (x = y = z = p x = p y = p z = 0) us- 
ing the computer algebra platform "Mathematica" (Wol- 
fram Research Inc.). The program allows us to compute 
the coefficients of the expansion using arbritrary precision 
arithmetics, while in Paper I the corresponding expansion 
was made with less accuracy. This change improves the 
credibility of our computation. After the Taylor expan- 
sion the Hamiltonian becomes H = X^=2 where H s is 
a homogeneous polynomial of degree s in x, y, z,p x ,p y ,p z . 
The second order term is 

Hi = \{pl+pl)~2xp y + {\ 



2" 

3\/3m 



— )x 



9^ 2 



-xy 



(3) 
(4) 



The last change of variables 
(x,y,z,p x ,p y ,p z ) -> (xi,X2,x 3 ,yi,y2,y 3 ) 
is performed in order to bring the quadratic part H2 
(Eq. ||) to the diagonal form 



Ho 



1 



(5) 



Since the term \{p 2 z + z 2 ) in Eq. (|^) has this form with 
lj 3 — 1 , the canonical transformation (Eq. ^) must bring 
to diagonal form the part of H 2 which corresponds to the 
planar case and depends only on x, y, p xi p y . This trans- 
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formation is the one performed in Paper I, but we include 

it to make the paper self-consistent. So we have 

(x,y,p x ,p y ) T = C {x 1 ,x 2 ,yi,y2) T , (6) 

where T denotes the transpose matrix and C is the matrix 

_i _i _i _i 
C ={e 1 m 1 2 ,e 2 m 2 2 ,f 1 m 1 2 J 2 m 2 2 ) (7) 

with 

+ 4V3a + 9 4a + 3%/3 4y^a + 9\ 



U= 0,2^- 



().■ 



Uj (8w| + 4y/3a + 9) Uj (4a + 3 ' 




27 

— +4a 2 
4 



(1 - 2^)3\/3 



We remark that the transformation of Eq. (|TJ) in- 
volves only the variables x, y, p x , p y since z, p z remain 
unchanged. For notational consistency we complete the 
transformation of Eq. (||) by putting x 3 — z and 2/3 = p z . 

All the above transformations were performed in order 
to bring the Hamiltonian to the form 

H(xi,x 2 , £3,1/1,2/2,2/3) = ^2H s (x 1 ,x2,X3,y 1 ,y2,y 3 ) ,(8) 

s>2 

where i/ s is a homogeneous polynomial of degree s in the 
hereafter called 'old variables' x\, X2, X3, 2/1, 2/2, J/3- The 
quadratic part H 2 which is given by Eq. (||) is the Hamil- 
tonian of a system of three harmonic oscillators with fre- 
quencies wi ~ 9.967575 • 10" 1 , u 2 ^ -8.046388 • 1(T 2 and 
— 1. The form of Eq. (11) is suitable for the direct ap- 
plication of the normal form theory, as it is described in 
detail by Giorgilli et al. fll989Q . We give a brief sketch of 
this procedure: we construct a generating function X, of 
the form 

A = ]TA S , (9) 

s>3 

where X s is a homogeneous polynomial of degree s, so 
that the corresponding canonical transformation brings 
the Hamiltonian to normal form 

(10) 



z = T z s , 



s>2 



where Z s is a homogeneous polynomial of degree s in the 
new 'normal variables' x' x , x' 2l x' 3 , y[, y' 2 , y' 3 . The term 
"normal form" means, that Z is a function of the quanti- 
ties I- = \{x^ + y t 2 ) with i = 1, 2, 3 , so that th e syst em 



is formally integrable (Birkhoff 1927; Gustavson 1966) 



The generating function X and the normal form Z are 
computed by solving the equation 

T X Z = H, (11) 
where Tx is an operator whose action is defined as follows 

T X Z = Y. F * > 



k>l 



where 



Fk — 2^ Zs,k- 



Z s ,o — Zs 



k 

Zs > k = Hk' Lx2 

n—1 



J s,k—n 



and Lx k Z s = [Xk,Z s ], with [ , ] denoting the Poisson 
bracket. The operator Tx is linear, invertible and pre- 
serves products and Poisson brackets (Giorgilli & Galgani 
1978 ). The computer program that solves Eq. ( |TT| ) and de- 



termines the generating function X and the normal form 
Z is described by Giorgilli fll979| ). 

Since the Hamiltonian in the old variables (Eq. ^) is an 
infinite series, in practice we stop the expansion at some 
order f and use the truncated Hamiltonian 

=H 2 + H 3 + --- + H f . (12) 
Then for any fixed integer r with 3 < r < f we solve 



Eq. (11) defining a truncated generating function X^ up 



to order r 

A« = X 3 + Xi + ■ ■ ■ + X r (13) 
and constructing the normal form Z^ up to this order 

z^ = z 2 + z 3 + --- + z r + y« , ' (14) 

where is a remainder, actually a power series starting 
with terms of degree r+1. So we have the equation 



T x^) H 



(r) _ 



z 



3 



z r 



normal form 



y{r) 



1 • ma null / 



(15) 



(r) 

The first term of the remainder is also computed, 

since it is needed for the estimation of the size of the 
effective stability region as we will explain in Sect. 3. 

In Paper I, where the planar problem (2D case) was 
studied, the power series of 4 variables were truncated at 
order f = 35. A function of 4 variables expanded up to 
order 35 requires 82,251 coefficients, while the process of 
constructing the normal form requires the computation 
of several functions with a total of 2,549,782 coefficients. 
In the spatial problem (3D case) we use expansions of 
functions of 6 variables up to order f = 30. This is a much 
harder task compared to the 2D case since a function of 
6 variables expanded up to order 30 requires 1,947,792 
coefficients and the program which calculates the normal 
form manipulates 55,929,459 coefficients. 

3. Estimation of the effective stability region 

3.1. Theoretical framework 

The transformed Hamiltonian Z^ admits three approxi- 
mate first integrals of the form 

I' J {x\y')= l -(x' J 2 +yf) , j = 1,2, 3. (16) 
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The variation of these quantities in time, is given by 



i' J = [Ij,Z^] = [r j ,Y^], j -1,2,3, (17) 
which is a power series starting with terms of degree r + 1. 

We introduce now suitable domains in the phase space, 
where we study the stability properties of the system and 
also a norm in these domains, in order to estimate the 
time variations of the three approximate integrals /j. 

For arbitrary fixed positive constants R±, i?2, R3 we 
consider a family of domains of the form 



A p r. = 



,'2 



y? 



< 



(18) 



where p is a positive parameter and x' , y' stand for x[, 
x' 2 , x' 3 and y[, y' 2 , y' 3 respectively. From the definition of 
the domains A p r it is evident that 

(x',y') G A p fl^/; < \ P 2 R) , j 



1,2,3 



(19) 

The norm H/Hpii of a homogeneous polynomial 
f(x',y') of degree s in the domain A p r does not exceed 
the quantity 

ll/IU< 

A 

2 s / 2 



IC 



jljtjaklkaka 



(20) 



J v>3 

1, 2, 3. We remark that for the 



E 

as shown in Paper I. Cj 1 j 2 j 3 k 1 k 2 k a are the complex co- 
efficients of f(x',y') when / is transformed in complex 
variables £, r\ via the transformation = (£j + irjj)/ \/2 

V'j = i(£j-iVj)/V2ior j 
above norm the elementary property 

||/IU = p1I/I|h 

holds. 

We assume that the initial point of an orbit lies in the 
domain A Po r for some positive value po, and we require 
that the orbit should be confined inside a domain A p r 
with p > po for a finite time interval. We shall refer to 
this time interval as the escape time r. Since I'j = dlj/dt, 
we get 



(21) 



dt > 



dl'j 



J = l,2,3, (22) 

su Pa p « 
where sup ApR 

Ij , or in other words the supremum norm of Jj , over the 
domain A p r. 



3=1,2,3 2 (r 



-M - 



1 



1 



(24) 



l)ll[^,n+Ullfl 
The above quantity depends on the order r up to which 
the normal form is constructed and on the radii po and 
p of the initial and final domain respectively. In Eq. (24) 
we keep the smallest value with respect to j (j = 1, 2, 3), 
because when the orbit is outside the disk with radius pRj 
in one of the planes (x^-, t/') it is also outside A p r. 

In order to have the minimum escape time as a func- 
tion of po we eliminate the dependence of T r (po,p) on p 
and r. It is evident that for a given order r the escape time 
goes to infinity as the radius of the outer domain grows. 
So by fixing p to be equal to Apoj with A > 1 we get 
T r ,x{po) = 

1 



1 - 



i=i*3 2(r-l)pTW,l&]ll* ' X ' 1 
By giving a fixed value to A we eliminate the depen- 
dence of the minimum escape time on the radius of the 
final domain. In particular, we put A = 1.2, which means 
that the radius of the final domain is 20% greater than 
the radius of the initial domain. 

The next step is to optimize the minimum escape time 
with respect to r. For A = 1.2 we compute T ri i.2(po) via 
Eq. ( |25| ) for r running from 3 to the maximum order f—1, 
for every value of po- We choose the optimal order r opt of 
the expansion as the one that gives the maximum value of 
the escape time. Thus we get the maximum escape time 
T as function of only the radius po of the initial domain: 
T(p ) = max r r , 1>2 (po) ■ (26) 

3<r<r 

The maximum order of the power expansions in the 
3D case is f = 30, which means that the normal form 
was computed up to order 29 and the first term of the 
remainder is of order 30. So it becomes clear that r opt < 
29. 



"31 



denotes the maximum absolute value of 3.2. General results 



Giorgilli ct al. (1989) proved that the power series of 
the remainder is absolutely convergent in a domain 
A p f{ provided p is small enough. So, assuming that p is 
smaller than half of the convergence radius of the remain- 
der we can use the approximate estimation 

(23) 



sup |Jj| 

(r) 

where V +1 



<2||[/',y; ( ; ) 1 || pJ , 2p 



(21) 



(r) 



r+1 



\R 



Y. 



(r) 



is the first term of the remainder. The term 

r + 1 and it 



r M is a homogeneous polynomial of degree 
is easily computed as a byproduct of the program that 
calculates the normal form. The validity of the above as- 
sumption is discussed in Sect. 4. 

By integrating both parts of Eq. ( p2[ ) and using 
Eq. (|23|), we estimate the minimum escape time as 

Tr(P0,P) = 



For a general discussion and for making the results com- 
parable to the ones found in Paper I for the 2D case, we 
put R± = i?2 = R3 = 1. In Fig. [I] we plot the logarithm of 
the maximum escape time (log T) as a function of the log- 
arithm of the radius of the initial domain (logpo)i m the 
spatial (3D) and the planar (2D) cases. In both cases the 
normal form is computed up to order 29. Taking an ini- 
tially small domain around the L4 point, all the orbits are 
confined inside a slightly larger domain (with 20% greater 
radius), for very long time intervals, while for large initial 
domains (large values of po), the escape time is small. 

A meaningful time interval for the system is the esti- 
mated age of the universe, which is in our time units 10 10 . 
This value is marked in Fig. [I] by a horizontal line and it 



corresponds to logpo 



■1.625, i.e. p ~ 2.371 • 10" 



the 3D case and to logp ~ -1.570, i.e. p ~ 2.692 • 10~ 2 
in the 2D case. In both cases the optimal order of the ex- 
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Table 1. Estimated stability region for some real asteroids. CN is the catalog number of the asteroids and Ri, R2, 
i?3 are the radii corresponding to the initial data. The radius po of the effective stability region is computed in the 
spatial (3D) case where the normal form was constructed up to order 29 and in the planar (2D) case where the normal 
form was constructed up to order 49. An asteroid is inside the stability region if po > 1. The optimal order of the 
expansion r opt is also reported in both cases. The table is sorted in decreasing order with respect to the value of po 
for the spatial case. 



3D case 2D case 



CN 


Ri 


R2 


R3 


Po 








Po 






fopt 


89211605 


0.033150 


0.019594 


0.013270 


1.0105 






29 


1.1722 






35 


2357 


0.042346 


0.028509 


0.049056 


6.3957 


10" 


1 


29 


8.7950 


10" 


1 


36 


88181612 


0.031302 


0.002101 


0.068011 


6.3294 


10" 


1 


28 


1.5364 






33 


41790004 


0.016517 


0.031063 


0.068685 


6.2154 


10" 


1 


29 


1.1687 






36 


5257 


0.031836 


0.042424 


0.031626 


6.1978 


10" 


1 


29 


8.0918 


10" 


1 


38 


4867 


0.233120 


0.736264 


0.468336 


4.0589 


10" 


2 


29 


5.1663 


10" 


2 


36 


1867 


0.216926 


0.758254 


0.466126 


3.9819 


10" 


2 


29 


5.0305 


10" 


2 


36 


88172500 


0.242971 


0.902082 


0.524644 


3.3747 


10" 


2 


29 


4.2302 


10" 


2 


36 


2363 


0.293736 


1.012524 


0.553082 


3.0026 


10" 


2 


29 


3.7667 


10" 


2 


36 


1208 


0.361925 


0.997579 


0.543939 


2.9917 


10" 


2 


29 


3.7807 


10" 


2 


36 




-5 I ' 1 1 1 

-2.5 -2 -1.5 -1 -0.5 

log po 

Fig. 1. The logarithm of the maximum escape time log T 
as a function of the logarithm of the radius of the initial 
domain logpo, in the spatial 3D case (solid line) and in 
the planar 2D case (dashed line) . In both cases the normal 
form has been computed up to order 29. The horizontal 
line marks the time that corresponds to the age of the 
universe. 



.4 I , , , , 1 

5 10 15 20 25 30 

fopt 

Fig. 2. The logarithm of the radius log po of the effective 
stability region which ensures stability for time equal to 
the age of the universe, as a function of the optimal order 
fopt of the expansion of the normal form, in the spatial 3D 
case (solid line) and in the planar 2D case (dashed line). 
In both cases the normal form has been computed up to 
order 29. 



pansion is r opt = 29. The best previously found estimation 
of the radius of the effective stability region, was obtained 
in Paper I in the 2D case, namely po ~ 2.911 • 10~ 2 , where 
the normal form was computed up to order 34. In that 
case the optimal order was r opt = 34. The radius of the 
effective stability region in the spatial case is about 12% 
smaller than the one computed for the planar case for 
r opt = 29 and about 19% smaller than the one found in 



Paper I for r opt — 34. So, the estimated stability region in 
the 3D case is a realistic one since it is comparable to the 
region found in Paper I. 

The estimated radii in the 3D and the 2D cases, for 
the same order of expansion of the normal form, are rela- 
tively close to each other since they are of the same order 
of magnitude, with the radius computed in the 3D case 
being always smaller, as seen in Fig. § Thus the Arnold 
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diffusion, which appears only in the 3D case, does not af- 
fect the size of the effective stability region significantly. 

We notice in Fig. || that up to order r opt — 14 the 
increment of the order r opt improves the estimation of the 
radius significantly both in the 3D and the 2D case. For 
orders greater than 15 the increment of the order leads to 
big increment of the computational effort needed for the 
construction of the normal form (mainly in the 3D case), 
but to relatively small improvements for the estimated 
radii. For instance, the radius found in the 2D case for 
fopt = 13 is almost equal to the one found for r op t = 14 
in the 3D case (logpo ~ —1.84) which was obtained by 
computing almost 20 times more terms for the various 
expansions, than in the 2D case. As r opt increases things 
become even more difficult. For r opt — 29 in the 3D case we 
find almost the same radius as in the 2D case for r opt = 25 
by computing almost 86 times more terms compared to the 
2D case. So, it is evident that pushing the computation of 
the normal form to higher orders becomes impractical for 
the 3D case. On the other hand this can be done more 
easily in the 2D case. 

According to Nekhoroshev theory, the series arising 
from the classical perturbation theory have an asymptotic 
character in the sense that one should not exceed an op- 
timal value for the order of expansion, which gives the 
best possible result. From Fig. || we see that this limit is 
greater than 30, both for the 3D and the 2D case. In Paper 
I, where the construction of the normal form was done up 
to order 34 for the 2D case, the limit r op t was not reached. 

By computing the normal form up to order 49 for the 
planar problem, so that the first term of the remainder 
is of order 50, we find the maximum value of the optimal 
order of the expansion of the normal form to be r op t = 38, 
and the corresponding value of logpo equal to —1.506, i.e. 
po ~ 3.119 • 10~ 2 . This value is about 7% greater than 
the one found in Paper I for r opt = 34. This is the best 
estimation of the radius of the stability region one can 
achieve using the particular theoretical framework, since 
expansions of the normal form to orders greater than 38 
do not improve the results. 

3.3. Application to real asteroids 

In order to apply the above results to the real solar sys- 
tem we examine if 98 real asteroids, which are located 
near the Lagrangian point L4, are inside the estimated 
effective stability region. Using the same catalog as in 
Paper I we extract the elements of the Trojan asteroids 
at the epoch, December 14, 1994, J.D.= 2449700.5. We 
also find the elements of Jupiter at the same epoch. As- 
suming that the orbit of Jupiter is circular we find the 
position of the L4 point and compute the coordinates 
{Qi, Q2, Q3, Pi, P%, P3) of the asteroids in the rotating he- 
liocentric system with the z axis orthogonal to the plane 
of Jupiter's orbit. By using the canonical transformations 
described in Sect. 2 we find the position of every aster- 



oid in the (xi, x%, X3, y±, j/2, IJ3) coordinates which diag- 
onalize the quadratic part (Eq. |3|) of the Hamiltonian. 
Then we define the radii R\, R2, R3 for every asteroid 
as Rj — ^Jx 2 + y 2 for j — 1,2,3. Using these values we 
determine the radius po of the effective stability region in 
the way described previously in this section. It is evident 
from the definition of R\ , R2 , R3 , that an asteroid is inside 
the stability region if po > 1. 

In the above procedure we use exactly the same aster- 
oids and their elements at the same epoch, as in Paper I. 
In this way our results can be compared directly to the 
ones obtained in Paper I where the effective stability of 
4 real asteroids was guaranteed. We apply this procedure 
both to the spatial case, where the normal form is com- 
puted up to order 29 and to the planar case, where the 
normal form is computed up to order 49. We remark that 
in order to find the radius po in the 2D case only the values 
of R\ and Ri are used, since we project the real asteroids 
on the plane of Jupiter's orbit. 

In Table [l] only the results for the 5 best and the 5 
worst cases of the radius p$ in the 3D case are presented. 
For every asteroid the corresponding results in the 2D case 
are also reported. In the 3D case one real asteroid is inside 
the stability region, while in the worst case the estimated 
value of the stability region's radius is smaller by a factor 
34. In all cases the optimal order of the normal form is 
the maximum possible, r opt — 29, which means that the 
results may be improved for higher orders. 

In the 2D case the results are slightly better than the 
ones achieved in Paper I, mainly because we performed 
the expansion of the normal form up to a higher or- 
der. This improvement does not change the results sig- 
nificantly. Four real asteroids (three reported in Table 1 
and one not shown in that table) are inside the planar 
stability region. In the worst case for the planar problem 
(asteroid 2363) a factor 27 is needed in order for the as- 
teroid to be safely inside the stability region. The optimal 
order for all asteroids is r opt < 38, although the expansion 
of the normal form was performed up to order 49. So the 
computation of the normal form to orders higher than 38 
does not improve the estimations in the 2D case. 

We remark that one would expect to find fewer as- 
teroids inside the stability region in the 3D case than in 
the 2D case, since the spatial stability region is projected 
on a plane in the 2D case. So points that are outside the 
spatial stability region may be projected inside the planar 
stability region. 

4. Discussion 

In Sect. 3.1 we assumed that the supremum of |7j | , which is 
an infinite series, can be approximated by the norm of the 
first term multiplied by 2 (Eq. ^3|). The practical reason 
for doing this approximation is that we cannot compute 
the whole infinite series, while on the other hand its first 
term can be obtained easily. The above assumption for 
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the norm is valid if the estimated value po ~ 2.371 • 10 -2 
in the spatial case, is smaller than the half of the conver- 
gence radius of the remainder Y^ r \ But estimating the 
convergence radius of the remainder is not possible since 
it requires the computation of higher orders of the series. 

One can attempt an estimation of the convergence ra- 
dius following Paper I. Since we have computed the ex- 
pansion of the generating function X (Eq. |9|) up to order 
30, we can evaluate its convergence radius by fitting the 
norms ||X s ||.r for 3 < s < 30 with a geometric sequence, 
i.e. we look for constants a and b such that 
||X s ||^<a6^ 3 . (27) 

In order to satisfy this condition it is sufficient to set 

i 

a =11X311* , 6 = 3 max o f|^V" • (28) 



3||R 



Giorgilli ct al. (1989) proved that if the generating func- 
tion X satisfies Eq. (|27j) then the coordinate transforma- 
tion x' — Txx, y' — Txy is absolutely convergent in a 
domain A p /j provided p is small enough. The same was 
proven for the inverse coordinate transformation and for 
the transformation of any other function such as the ac- 
tions Ij (Eq. |TJ). So, we fit with a geometric sequence 
the norms of the coordinate transformations Txxi, TxX2, 
TxX3, Txyi, TxTj2, Txyz and of the transformations of 
the approximate integrals to the old variables I(, 
T^I 2 , T^I's- The worst case (T^ 1 ^) gives b ~ 18.665, 
which corresponds to a convergence radius p ~ 5.358T0 -2 . 
Considering this value as a good indicator of the true con- 
vergence radius, we conclude that the estimated radius of 
the stability region is smaller than the convergence radius 
of the series, by a factor 2.2, so it is safely inside the con- 
vergence domain and the approximation used in Eq. ( p3f ) 
holds. 

The theoretical framework we used proved to be very 
efficient, since we are able to guarantee the effective stabil- 
ity of one and four real asteroids in the spatial and planar 
problem, respectively. On the other hand we have reached 
the limits of its effectiveness, since the above results can- 
not be improved significantly. 

Computing the expansion of the normal form to higher 
orders in the 3D case would slightly improve the esti- 
mations, as we can see from the respective estimations 
performed in the 2D case, where the optimal order was 
reached for r opt = 38. 

The estimation of the minimum escape time T r {p ,p) 
(Eq. ^) was improved compared to the one used in Pa- 
per I. This improvement was done by taking into account 
the dependence of sup ApR on the radius p of the final 
region (Eq. |23| ), while in Paper I this quantity had been 
overestimated by considering sup ApR constant inside 
the domain A p # and equal to its value at the edge of 
the domain. The above change forced us to fix the ra- 
dius p of the final domain. In our calculations we used 
A = p/pa = 1.2 but this value does not influence strongly 
the results. This can be easily understood since, the quan- 



tity [1 - (l/A*"" 1 )] in Eq. (§|) is very close to 1 for A > 1.1 
and r = 29 (which is the order of expansion up to which 
we compute the normal form). The improved estimation 
of the escape time does not change the results drasti- 
cally. For instance, in the 2D case for Ri = R2 = 1 with 
the normal form being computed up to order 34, we got 
po ~ 3.005 -10~ 2 , which is about 3% greater than the value 
po ~ 2.911 ■ 10~ 2 obtained in Paper I for the same case. 

The norms YJ.+x]H-R for 3 = !>2,3 in Eq. (§5|) 
were estimated using Eq. (pc|). Since this estimation is 
purely analytic it is certainly pessimistic. For this reason 
we made numerical calculations of the norms in order to 
determine whether the overestimation of the norms influ- 
ence strongly the estimation of the radius of the stability 
region. We used two Fortran codes. Algorithm GLOBAL 



(Bocnder et al. 1982), which had the best performance in 
all cases, was used to obtain the norms, and algorithm 
SIGMA (TOMS 667) (Aluffi-Pentini et al. [l988a| , §) was 
used for verification purposes. GLOBAL is a stochastic 
algorithm for finding the maximum of a real-valued func- 
tion. In stochastic methods for optimization, the proba- 
bility of finding the global maximum approaches unity as 
the sample size of the random initial values increases. This 
algorithm utilizes a combination of sampling, clustering, 
and local search; it terminates with a range of confidence 
intervals on the value of the global maximum. SIGMA 
is a global optimization algorithm, which implements a 
method founded on the numerical solution of a Cauchy 
problem for a stochastic differential equation inspired by 
statistical mechanics. A global maximum of the function is 
sought by monitoring the values of the function along tra- 
jectories generated by a suitable discretization of a first- 
order stochastic differential equation. 

In all cases the numerically found maximum of 
|| [I'j , 3^.+]] || r, j — 1,2,3 was located at the edge of the do- 
main Aji. So by limiting our computation on the edge we 
were able to evaluate the maximum with even greater ac- 
curacy. The analytically estimated maximum was greater 
than the one computed numerically by a factor 28 in the 
worst case. Using the numerically calculated norms in 
Eq. ( p6| ) we found p ~ 2.667 ■ 10~ 2 . This value is only 
1.12 times greater than the one computed by using the 
analytically estimated norms. So the improvement of the 
size of the stability region is negligible compared to the 
factor 34 needed for all the real asteroids to be inside the 
stability region. Thus the estimations based on the ana- 
lytically found norms are reliable. 

5. Summary 

We studied the spatial and the planar circular restricted 
problem of three bodies in the spirit of Nekhoroshev the- 
ory, considering in particular the problem of the stability 
of the Trojan asteroids around the Lagrangian point L4 
in the Sun- Jupiter-Asteroid model. By constructing the 
normal form with terms up to order 29 in the spatial case 
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and up to order 49 in the planar case and using analytical 
methods and numerical approximations, we made realistic 
estimations of the size of the effective stability region. Our 
results can be summarized as follows: 



1. The estimated size of the effective stability region in 
the spatial case is big enough to include 1 real aster- 
oid. The region where the most remote asteroid is lo- 
cated (out of the 98 real asteroids we checked), is larger 
by a factor 34 compared to the estimated stability re- 
gion. This result imp roves significantly ol der es timates 
(Giorgilli et al. |l989| ; Celletti & Giorgilli |l99l| ) where 
no real asteroid was inside the stability region and a 
factor 3,000 was needed for the most remote asteroid 
to be inside the stability region. 

2. The radii of the effective stability region in the general 
spatial and planer cases are close to each other for 
the same order of expansion of the normal form, with 
the radius computed for the spatial case being always 
slightly smaller. Thus, Arnold diffusion does not affect 
the size of the effective stability region significantly. 

3. By computing the normal form in the planar case, to 
higher order than in Paper I, the optimal order of the 
expansion of the normal form was actually reached and 
it was found to be r opt = 38. Also the estimation of the 
stability region's size was slightly improved compared 
to Paper I. Four real asteroids are found to be inside 
the stability region, as it was found also in Paper I. 

4. Several improvements in the computational procedure, 
compared to previous works, were introduced in the 
present paper, namely: a) the computation of the nor- 
mal form to higher orders; b) a more accurate estima- 
tion of the minimum escape time r r (po,p) (Eq. 



c) the numerical estimation of the norms ||[Jj, 5^.+\]||.r 
for j — 1, 2, 3, which also verified the reliability of the 
analytically calculated norms via Eq. (^o|) and d) the 
use of an arbitrary precision computer algebra system 
for the expansion of the Hamiltonian. 
5. The theoretical framework we used reached the limits 
of its effectiveness by providing the best possible re- 
sults in the 2D case, and comparable results in the 3D 
case. In order to have a non negligible improvement of 
the estimated size of the effective stability region, one 
has to choose better coordinates than the cylindrical 
ones used in this paper, in the sense that these coordi- 
nates should be more adapted to describe the banana- 
shaped region of the actual stability region around L4. 
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